Simulations of a single membrane between two walls 
using a Monte Carlo method 
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Quantitative theory of interbilayer interactions is essential to interpret x-ray scattering data and 
to elucidate these interactions for biologically relevant systems. For this purpose Monte Carlo 
simulations have been performed to obtain pressure P and positional fluctuations a. A new method, 
called Fourier Monte-Carlo (FMC), that is based on a Fourier representation of the displacement 
held, is developed and its superiority over the standard method is demonstrated. The FMC method 
is applied to simulating a single membrane between two hard walls, which models a stack of lipid 
bilayer membranes with non-harmonic interactions. Finite size scaling is demonstrated and used 
to obtain accurate values for P and a in the limit of a large continuous membrane. The results 
are compared with perturbation theory approximations, and numerical differences are found in the 
non-harmonic case. Therefore, the FMC method, rather than the approximations, should be used 
for establishing the connection between model potentials and observable quantities, as well as for 
pure modeling purposes. 



I. INTRODUCTION 

Recent research on lipid bilayers [Q has contributed 
to the important biological physics goal of understand- 
ing and quantifying the interactions between membranes 
by providing high resolution x-ray scattering data. From 
these data the magnitude of fluctuations in the water 
spacing between membranes in multilamellar stacks is 
obtained. This enables extraction of the functional form 
of the fluctuational forces, originally proposed by Hel- 
frich (|] for the case of hard confinement. For systems 
with large water spacings, the Hclfrich theory has been 
experimentally confirmed j|. For lecithin lipid bilayers, 
however, the water spacing is limited to 20A or less. For 
this important biological model system, our data show 
that a theory of soft confinement with a different func- 
tional form is necessary; this is not surprising because 
interbilayer interactions consist of more than hard-wall, 
i.e., steric interactions. 

The theory of soft confinement is even more difficult 
than the original Hclfrich theory of hard confinement. 
Progress has been made by modeling the stack of inter- 
acting flexible membranes by just one flexible membrane 
between two rigid walls ||||. Even with this simplifi- 
cation, however, the theory involves an uncontrolled ap- 
proximation using first order perturbation theory and a 
self-consistency condition in order that the interbilayer 
interaction may be approximated by a harmonic poten- 
tial ||. We have obtained inconsistent results when ap- 
plying this theory to our data (unpublished). Possible 
reasons are (i) the theory is quantitatively inaccurate or 
(ii) the single membrane model is too simple. The im- 
mediate motivation for this paper is to test possibility 

In order to obtain accurate results for a system with 
realistic non-harmonic potentials, we use Monte-Carlo 
(MC) simulations. The particular MC method developed 



in this paper will be called the FMC method because it 
uses the Fourier representation for the displacement of 
the membrane rather than the customary pointwise rep- 
resentation, which will be called the PMC method. The 
main advantage of the FMC method is that the opti- 
mal step sizes do not decrease as more and more ampli- 
tudes are considered. In contrast, in PMC simulations, 
the optimal step sizes decrease as the inverse of the den- 
sity of points in one dimension, because the bending en- 
ergy becomes large when single particle excursions make 
the membrane rough. Because of this, relatively large 
moves of the whole membrane are possible with the FMC 
method, but not the PMC method. This produces rapid 
sampling of the whole accessible phase space, while re- 
specting the membrane's smoothness. The resulting time 
series have moderate auto-correlation times Q that do 
not increase substantially as the membrane gets larger 
and/or more amplitudes are taken into account. Even 
though each Monte Carlo step takes longer, FMC still 
outperforms PMC by a wide margin. It then becomes 
possible to carry out substantial simulations on a stan- 
dalone workstation rather than a supercomputer JtJ and 
to obtain accurate results for a single membrane subject 
to realistic potentials with walls, and even for a stack of 
such membranes (to be described in a future paper) ^ . 

Section [l| defines the membrane model and thephys- 
ical quantities simulated in the paper. Section HI de- 



scribes the FMC method and also gives some important 
details that are used to speed up the code. In Section 



[V the method is tested on an exactly solvable model, 



namely, one that has only harmonic interactions with 
the walls. This test also allows examination of the sys- 
tem properties and the convergence of FMC results for 
an infinitely large, continuous membrane. In section ^ 
the FMC method is applied to a single membrane with 
realistic, non-harmonic interactions with the walls. Sec- 
tion VI makes a detailed comparison of the FMC method 
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and the standard PMC method. This section shows that 
the FMC method not only converges faster to average 
values for continuous membranes, but also gives smaller 
stochastic errors. Finally, section VIl compares simula- 
tion results with those obtained using the analytic first- 
order theory of Podgornik and Parsegian || and from 
experiment (T|. 



II. SINGLE MEMBRANE SYSTEM 
L 
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FIG. 1. A fluctuating single membrane, constrained be- 
tween two hard walls. 



At the atomic scale a lipid membrane is composed of 
complex lipid molecules and many simulations are per- 
formed at this scale p[-pT]|. However, for modeling the 
structure factor for low angle x-ray scattering (in con- 
trast to modeling the form factor), it is customary and 
appropriate [O-lLq] to model the membrane as an in- 
finitely thin flexible sheet as shown in Fig.|l|. 

The membrane undulates with instantaneous fluctua- 
tions in the z-direction, given by u(x,y), subject to pe- 
riodic boundary conditions. The model energy W is a 
sum of bending energy with a bending modulus K c and 
an energy of interaction with the walls, 



W 



2 



(Ait) dx dy + / w a (u)dx dy 



(1) 



Since each wall is a surrogate for a neighboring membrane 
in a stack, and since it is desired to obtain physical prop- 
erties per membrane, the interaction potential is given by 
the average w a (u) = (V(a + u) + V(a — u))/2 of the inter- 
actions V with each wall and the corresponding volume 
of the system per membrane is then aL 2 . For a separa- 
tion z between a wall and the membrane the interaction 
potential will be based on the standard form 



V(z) = A\e- Z / X 



H 



127TZ 2 



(2) 



where the first term on the right hand side is a repulsive 
hydration potential || and the last term is an approxi- 
mate, attractive van der Waals potential. The divergence 



in the van der Waals potential as z— >0 in Eq.(||) is quite 
artificial; physically, it is masked by stronger steric repul- 
sions at small z p6[ . This is corrected in this paper by 
including only a finite number of terms m ma x in a power 
series expansion of 1/z 2 about u = 0. It is shown later 
that a wide range of m max give nearly the same result, so 
m-max is not a critical parameter and power series suffice 
to represent the van der Waals potential satisfactorily for 
the most probable values of z but avoid including artifi- 
cial traps near the walls. Other forms besides Eq.(||) can 
be treated as well. 

The first important quantity, obtained directly from 
the simulation, is the mean square fluctuation a 2 in the 
water spacing. In Fig.[j], a 2 = u 2 {x, y), where the average 
is over both space and time. The second physical quan- 
tity is the pressure P that must be exerted on the walls 
to maintain the average water spacing a. The pressure is 
a sum of two components: P±, caused by collisions and 
equal to a temporal average of a delta-function-like in- 
stantaneous pressure, and P2, which is due to non-contact 
interactions with the walls, and that varies smoothly with 
time and position. A virial theorem argument can be 
used to compute P\ . The general result is 



P = 



N 2 k B T-2U 
2aL 2 



2aL 2 



u——dx dy 

ou 



dw(u, a) 
da 



(3) 



where Pi is the term in square brackets. The relative im- 
portance of Pi and P2 depends on the potential. If the 
potential is completely steric (hard wall), then P2 = 0. 
However, we have found that for the more realistic poten- 
tials considered in this paper P± is very small compared 
to P2 because there are very few hard collisions. 



III. FOURIER MONTE CARLO (FMC) METHOD 

The membrane displacement u(x, y) is represented by 
its Fourier amplitudes u(Q), 

where Q — (2irm/L, 2irn/L), N is the total number of 
modes in each dimension and — N/2 + l<m, n<N/2. Re- 
ality of the displacement u(x, y) is guaranteed by requir- 
ing u(—Q) — u*(Q). Also, note that u(Q = 0)^0 allows 
the center of gravity to fluctuate away from the midplane 
between the walls. 

Using the standard Metropolis algorithm, the simu- 
lation attempts to vary one Fourier amplitude, picked 
randomly, at a time. The initial step sizes, which de- 
pend upon Q, are determined using a simplified form 
of the analytic theory ||. After a certain number of 
Monte Carlo steps (MCS), step sizes are adjusted using 
Dynamically Optimized Monte-Carlo(DOMC) @. Step 
size optimization results in an acceptance-rejection ra- 
tio of about 1/2, thereby minimizing the autocorrelation 
time t. In practice, because the initial values are already 
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based on a reasonably good approximation, DOMC ad- 
justment does not significantly improve the efficiency. 

The change in bending energy in Eq.(|]) after attempt- 
ing a step in u(Q) is K c L 2 Q A /2 times the change in 
|u(Q| 2 , which requires little time to compute. In contrast, 
calculating the change in the interaction energy with the 
walls requires a real space representation of u(x, y) . How- 
ever, it is not necessary to use a fast Fourier transform 
(FFT) routine because the linearity of the Fourier trans- 
form requires only recomputing one Fourier term in order 
to update u(x, y). The time this takes is only 0{N 2 ) com- 
pared to 0(N 2 In N) for a standard FFT routine. Incre- 
mental addition errors are negligible for the longest runs 
when double precision is used; alternatively, one could 
perform FFT at long intervals to control such an error. 
The natural choice is made to approximate the interac- 
tion integral over the membrane by a sum over a set of 
equally spaced points (Li/N,Lj/N), with 0<i,j < N. 



IV. HARMONIC INTERACTIONS AND 
FINITE-SIZE SCALING 

To test the simulation code and investigate conver- 
gence to an infinite, continuous membrane, it is useful 
to consider a harmonic interaction energy. It is also use- 
ful to relate the parameters in the harmonic potential 
to those in Eq.(|l|) by expanding w a (u) to second order 
about u = 0, 



AXexp(-afX) 1 + 



H 



YI-kcl 1 



X- 



(4) 



so that the realistic Eq. ([j]) then takes the completely har- 
monic form 



^ = T J(V 2 u(r)) 2 d?r + 
B(a) 



u z (r)d 2 r + w (a)L z 



(5) 



where B = {A/X)er a / X - H/(2na 4 ) and w (a) = 
AXer aj x - H j (12-7ra 2 ) . The exact solution (valid for finite 
L and N/ L) for this harmonic model is 



r,2 



I 



L 2 ~~1 K c ( q 2 x +q 2 y ) 2 + B- 



(6) 



is written for the general case of realistic potentials and 
can then be applied when <7<§CA. As an example, con- 
sider a membrane with parameters N = 4, L = 700A 
and a non- harmonic potential with A = 1, H = 100 
(ttw = 2), A = 10A, K c = 1, T = 323K, a = 20 A, 
where Jl8| gives the units for A, H and K c used in this 
paper. The simulation gives a = 0.3394±0.0004A and 
P = 1.2877-10 7 ±200er3/cm 3 . In this case, a exact = 
0.33954A, and 



P = Ae~ a/X 



a 

2X 2 



H 



4 



1.28774-10 7 er ff /cm 3 , 



(8) 



again showing that simulation results are precise. 

The second usage of Eqs. (||) and (Q) is to obtain a 
and P as functions of N and L through the finite sums 
over Q. Simulations are always done with a finite num- 
ber of Fourier amplitudes and a finite-sized membrane. 
However, real membranes are continuous and the rele- 
vant size may be larger than 1/ito. So it is important 
to see how the results for finite systems can be used to 
obtain quantities for dense (N— >oo) and large (L— *oo, 
N/L = const) systems. Eqs. (0) and (0) can be used to 
compute a(N,L) and P{N,L) numerically to examine 
the asymptotic behavior of these functions. The result of 
such analysis is an asymptotic relation 



Ci 



-C 2 



1 

L2 : 



(9) 



and 



where typically d - lO^A" 1 and C 2 - 10 3 A 3 . The 
variability caused by the C 2 term is very small; typically 
about 0.2% when L>700A. However, the Ci term causes 
a for a finite membrane to vary with N as much as 20%. 



V. OBTAINING RESULTS FOR REALISTIC 
INTERACTION POTENTIALS 

Tabic | shows results for two selected non-harmonic 
potentials and a variety of sizes. One may first note that 
the autocorrelation times r a i and rp are nearly constant 
with system size. Next, convergence with increasing N 
and constant L is shown in Fig.|2j when the vdW interac- 
tion is absent. This behavior is similar to that of a har- 
monic interaction. The limiting values can be estimated 
by fitting the curve y = yoc + C 2 /N 2 + C 3 /N 3 . The fits, 
shown as solid lines on Fig.||, lead to cr^ = 4.394±0.004A 
and Poo = 202400±700er 5 /cm 3 . 



P = Ae~ a/X 



a 

2X 2 



(7) 



Equations (Q) and (^) are useful in two ways. First, the 
harmonic approximation given by Eq.|] is good if cr<cA. 
That provides a test of the correctness of the code, which 



3 



TABLE I. Representative simulation results for two interactions. 



N 


L[A] 




MCS,10 3 




TP 






A = 1, H = 0, K c = 1 [L8|, A = 


1.8A, T = 323if, a = 20A 







4 


700 


4.0774±0.0018 


123010±170 


500 


1.59 


1.35 


6 


700 


4.2767±0.0034 


156100±400 


100 


1.44 


1.18 


8 


700 


4.3376±0.0028 


173700±400 


100 


1.19 


0.96 


8 


700 


4.3366±0.0013 


173470±170 


500 


1.21 


0.98 


12 


700 


4.359±0.008 


187000±1300 


10 


1.16 


0.97 


16 


700 


4.3792±0.0034 


193800±600 


50 


1.08 


0.88 


24 


700 


4.3864±0.0024 


197920±430 


30 


0.946 


0.768 


32 


700 


4.399±0.011 


201500±1900 


6260 


1.43 


1.41 


32 


700 


4.3976±0.0030 


200600±500 


20000 


0.955 


0.741 






A — 1, H = 3, t7i max 


= 4, K c =0.1, A = 


1AA, T = 323K , a 


= 17A 




4 


350 


6.0902±0.0027 


28000±900 


500 


2.46 


1.03 


6 


525 


6.1097±0.0029 


34400±900 


200 


2.74 


0.96 


8 


700 


6.1225±0.003 


38500±1000 


100 


2.7 


0.97 


12 


1050 


6.128±0.005 


40800±1500 


20 


2.73 


1.05 


16 


1400 


6.1270±0.0026 


40000±600 


30 


2.35 


0.86 


32 


2800 


6.136±0.003 


42000±600 


6 


2.65 


0.89 



4 




0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 
1/N 2 

FIG. 2. a and P vs. 1/N 2 for A = 1, H = 0, A = 1.8A, 
a = 20l, A" c = 1, T = 323K and L = 700A. 
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I I -1 

0.00 0.01 0.02 0.03 

1/N2 

FIG. 3. a and P(1/N 2 ,L = const = 700A) for A = 1, 
H = 3, m maa; = 4, A = 1.4, a = 17 and K c = 0.1. The lines 
are drawn to guide the eye. 

Unfortunately, one does not obtain the same asymp- 
totic behavior as in FigJ| when the attractive force is 
large enough that the total potential has a maximum 
rather than a minimum when in the middle of the space 
between the walls. For instance, when Hj^O, a first de- 
creases with N, although later it gradually levels off and 
appears to have a minimum. It is interesting that, while 
a may change in an unexpected way as N increases, for 
the interaction considered, the pressure is still a smooth 
quasi-linear function of 1/N 2 (N— >oo), as shown in Fig. 
y, and its limiting value as N— >oo can still be estimated 
by extrapolation. Despite these variations in convergence 
behavior, the associated changes in a become very small 
and are certainly less than the desired accuracy of 1-2%, 
so we suggest that it is sufficient to increase N only to 
the point where further increases result in changes in a 
and P that are less than the target precision. 



The other variable that is potentially significant is the 
size of the membrane. Any physical quantity may depend 
on how large the membrane is, attaining a certain lim- 
iting value as L— >oo. By increasing L while keeping the 
"density" N/L — const, the membrane size is determined 
for which a and P approach their limiting values suffi- 
ciently closely. As in the case of harmonic interaction, 
the changes in these quantities are relatively small as L 
is increased. Indeed, when there is no attractive force, 
the changes are so small that they cannot be resolved 
reliably even when the estimated statistical errors are of 
order of 3TCP 3 A. When the interaction is smaller, the 
trends become more pronounced and similar to those seen 
for the harmonic potential. An example is given in Fig.|] 
which shows that for a moderate sized membrane the re- 
sults approach smoothly and closely those for an infinite 
mejmbrane(-L— >oo). For L = 700 A the difference between 
the estimated limiting value of a and the observed one at 
700 A is less than 0.5%, while for the pressure the same 
difference is less than 5% which is about the same as the 
experimental uncertainty in P. 
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2.5 1 ' 1 ' 1 ' 1 ' 1 1 

0.0 2.0x10-6 4.0x10-6 6.0x10-6 8.0x10-6 

1 / L2[A] 

FIG. 4. cr and P vs. 1/L 2 with N/L = 8/700A for A = 1, 
H = 3, K c = 0.1 @, m max = 4, A = 1.4 and a = 17. 

In summary, of the two factors that could affect con- 
vergence of simulation results, i.e. N and L, N is most 
important. L is therefore fixed, typically at 700A N is 
increased until the changes in quantities of interest are 
less than the target precision. We then fit a simple func- 
tion such as y = j/oo + C2/N 2 + c s /N 3 to the sequence of 
finite N results to estimate y^. 
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VI. COMPARISON OF FMC AND STANDARD 
PMC METHODS 

A. Basics of the PMC Simulation Method 

The standard way to simulate membranes [Q will be 
called the pointwise MC (PMC) method in which the 
potential of the system is given in discretized form 

ij nn 

^E^K). ( 10 ) 

ij 

where ^2 nn u is the sum of displacements of nearest 
neighbors of site For a harmonic potential, w(u) = 

wq + Bu 2 /2, and for periodic boundary conditions the 
exact solution for the mean square displacement is 

a 2 = kBT^ {B + AKc ^ {cos{Qx L ) + 
Q 

cos(Q y ^)-2fy\ (11) 

where Q x>y = 2im/L, -y + l < n <f ■ As with the FMC 
method, such an exact solution is useful in checking cor- 
rectness of the simulation code. 

The standard Metropolis algorithm is used, moving 
one point at a time in the PMC method. To start the 
simulation, an effective B is estimated using perturbation 
theory ||. It is then used in a formula that gives the 
mean-square fluctuation of a point (assuming harmonic 
potential) about its equilibrium position, determined by 
its environment: 

aiocal = \j BL 2 /N 2 + 2QK C N 2 /I? (12) 

Eq.(|l2"|) gives the initial step size. After a certain number 
of steps, DOMC [[HJ is used to compute the optimal step 
size, which is used thereafter. Some results using the 
PMC method are presented in Table ||. 



TABLE II. Real space simulations of membranes with dif- 
ferent density of points, constrained by a harmonic potential 
with B = 8. 303- 10 11 erg /cm 4 obtained from A = 1, H = 0, 
K c = 1 (Hi, A = 1.8 A, a = 20A. T = 323K, L = 700A. 
Simulation lengths are measured in 10 6 MCS. 



N 


a [A] 


MCS 


MCSo.i%* 




4 


8.390±0.005 


1 


0.41 


4.36 


6 


8.481±0.008 


1 


0.98 


13.8 


8 


8.332±0.031 


0.2 


2.77 


41.9 


8 


8.347±0.032 


0.2 


2.94 


42.3 


8 


8.305±0.010 


2 


2.73 


39 


12 


8.073±0.016 


4 


14.9 


203 


12 


8.070±0.015 


4 


14.6 


198 


16 


8.00±0.06 


1 


66 


782 


16 


8.07±0.06 


1 


59 


709 



* A simulation of approximately such length would have 
to be done to attain 0.1% accuracy for a. 



G 



B. Comparison of the FMC and PMC methods 

The time required to obtain a target error is one of the 
issues determining the viability of any simulation tech- 
nique. It is impacted by two separate factors: the rela- 
tive magnitude of random errors, and the speed at which 
various quantities, obtained for a finite system, converge 
to their values for the continuous infinite system. These 
factors are now considered in detail, to demonstrate the 
improvements of the FMC method. 
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FIG. 5. Variation with N of the simulation length 
MCSq,x%, required for 0.1% precision of a, for a PMC sim- 
ulation of a harmonic potential with A = 1, H = 0, K c — 1, 
A = 1.8A, a = 20 A, T = 323^ and L = 700A and for an 
FMC simulation for a realistic model potential with the same 
parameters. 



The random errors in estimated averages depend on 
the autocorrelation times of generated time series. These 
times are an indication of how "natural" the chosen ba- 
sis is for the simulated system. In the case of harmonic 
interactions, the variables used by FMC are exactly inde- 
pendent and therefore it is possible to vary each of them 
separately over its whole range. Although they do be- 
come correlated for realistic interactions, one would still 
hope that their dependencies are not great, and so they 
still represent a good basis. For PMC simulations, how- 
ever, the motion of any point is constrained by its envi- 
ronment, so one would expect the quality of time series to 
deteriorate as the "density" of the membrane and the im- 
portance of the local environment increase. These asser- 
tions are supported by Tables || and ||, which show that 
for FMC the autocorrelation times remain roughly con- 
stant with increasing A, whereas for PMC r a increases as 
./V 4 . A related question is how the simulation length (in 
MCS) required to obtain a certain accuracy (chosen to be 
0.1%) varies with N. A straight line fit to ln(MCS 0A %) 
vs. In A dependence for PMC has a slope of approxi- 
mately 4 (Fig||) . Therefore, the amount of time required 
to obtain a with the same precision grows as A 6 for PMC 
method. A somewhat surprising result is that the length 
required to achieve a given error estimate with FMC de- 
creases with N (Fig]5j). The precise law governing this 
decrease is unclear because of the difficulty of estimat- 



ing autocorrelation times; one guess, supported by the 
four points in the middle (N = 8 through 24) is that 
the length decreases as 1/%/A"; however, the hypothesis 
of the length staying asymptotically constant cannot be 
ruled out either. Because 1 MCS (for FMC) takes the 
amount of time 0(N 4 ), the computational complexity of 
the process generated by a Fourier-space simulation is 
only A 3 - 5 or A 4 , assuming that the same error estimate 
is achieved. This is a significant improvement over the 
A 6 law for the real-space simulations. 

The second factor favoring FMC concerns how closely 
the bending energy is approximated by the discrete ap- 
proximation in Eq.(10). This can be evaluated by the 
exact result for a for a harmonic model. Fig.^ shows 
that one requires larger A to obtain the same precision 
with the discrete approximation to the bending energy 
required by the PMC method in Eq.(|l0|) than for the 
true continuum model that can be treated naturally by 
the FMC method. 



•< 




N 

FIG. 6. Exactly computed a(N,L = 700A) for 
Fourier-space (Eq.(g)) and real-space (Eqs.(flo|) and ([fl|)) 
models of a harmonic potential with B = 8.303-W 11 erg /cm 4 . 
The other parameters are K c = 1, T = and L = 700 A. 



A specific example illustrates the preceding principles 
and also gives some typical computer times for these sim- 
ulations. The example is the harmonic model with pa- 
rameters given in Figj^. For the PMC simulation, A = 46 
was chosen so that <J exact (46, L — 700A) = 7.7898 was 
within 0.5% of its value 7.7478A for a continuous mem- 
brane. A simulation of 800,000 MCS took 9.5 hours on an 
SGI workstation with MIPS R5000 1.0 CPU and 128 Mb 
of RAM, running IRIX 6.2 and resulted in a = 7.33±0.19. 
So, 9.5 hours were insufficient to obtain a with 0.5% 
accuracy, and about 9.5-(0.19/(0.005-7.75)) 2 «229 hours 
would be required to achieve that precision. Turning 
to FMC, for A = 16 the exact a = 7.711ll. A run 
of 10,000 MCS yielded a = 7.7184±0.0165 and required 
only 240 seconds on the same computer as the PMC sim- 
ulation. One may also compare the time it takes to ob- 
tain the same estimates of random errors for the same A 
for the two methods. To do this, N = 16 and a target 
error of about 1% were chosen for the same interaction 
as before. A PMC simulation for 300000 MCS took 1174 
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seconds on an SGI workstation with a similar configura- 
tion to the one used in the previous test and resulted in 
a = 8.032±0.082A (t e = 14.7, r CT = = 441), a slightly big- 
ger error than desired. In contrast, an FMC simulation 
(also with N = 16) for 2000 MCS took only 63 seconds on 
the same computer, and resulted in a = 7.674±0.070A 
(te = 2.19, t ct 2 = 1.44), the random error in a now being 
slightly better than the target. So, in addition to a much 
faster convergence of the expected value to one for a con- 
tinuous membrane, the FMC method is also the faster 
one to obtain a given estimate of stochastic errors. 



VII. RESULTS AND IMPLICATIONS 

A. Distribution of the membrane displacements 

The functional form of the probability density function 
(pdf) is a central assumption in the perturbation theory 
0. Also, the behavior of the pdf near the walls is sig- 
nificant in discussing the formal divergence of the van 
der Waals potential and the importance of the hard wall 
collision pressure Pi. If the pdf does not decay to zero 
sufficiently quickly near the walls, then the value of m max 
used in the power series expansion would be a sensitive 
parameter and one would expect many hard collisions 
with the walls. The inset to Fig.0 shows that the pdf de- 
cays to zero near the walls in much the way that is pos- 
tulated by theory ^ . This is consistent with our results 
that Pi is small and m max is an insensitive parameter. 
This latter point is explicitly illustrated in Fig.|| which 
shows that the results for a plateau for 6 < m max < 40; 
a similar plateau occurs for P. Finally, Fig.^ shows that, 
away from the walls, the pdf is noticeably different from 
the theoretically assumed pdf || and it is generally dif- 
ferent from a Gaussian. 
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Gaussian 
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FIG. 7. Membrane pdf for a realistic constraining poten- 
tial. A = 0.2, H = 0.5, A = 1.3A, m max = 3, T = 3237^, 
K c = 0.1, a = 22A, N = 32 and L = 700A. Also shown are 
the Gaussian pdf, corresponding to a = 8.0196A and the ap- 
proximate pdf for the case of pure steric constraint proposed 
in §(Eq.(20)). 




m 

max 

FIG. 8. The relationship between the number of terms in 
the expansion approximating van der Waals potential and a, 
for the parameter set a = 1, H = 6, A = 1.8, K c — 0.2, 
T = 323, a = 13, L = 700. The line is drawn to guide the 
eye. 



B. P and a 



Simulation results: 
■ H=4 o H=0 

P 2 theory 

H=4 H=0 




14 16 18 20 22 

a [A] 

FIG. 9. o(a) and \nP(a), obtained from a simulation for 
A = 1, H = 4, A = 1.8, K c = 0.2 and also for H=0 (all other 
parameters being the same) and corresponding results from 
the perturbation theory 



For any kind of interaction, the main results to com- 
pare to experiment are the relationships between InP 
and a, and a and a. Figure ^ shows InP and a for sev- 
eral values of a. Two interaction types are considered: 
A = 1, H = 4, A = 1.8, K c = 0.2 and the same set with 
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H = 0. These figures also show the results obtained from 
the first-order perturbation theory || • The largest differ- 
ences with the simulations occur at larger a and when H 
is non-zero. In particular, the theory under-predicts the 
value of a at P = when no osmotic pressure is applied. 
Overall, however, the theory predicts trends quite well. 



C. Comparison to Experiment 

Recently, it has been proposed that the pressure due to 
fluctuations, Pfi, can be obtained from x-ray line shape 
data jjj. The derivation involves the use of harmonic 
Caille theory (l2Jl5|], which yields 



fi 



4ksT 

7T 8 



1 da~ 2 
K r da 



where a is obtained from 



t] 1 D 2 /tt 2 



(13) 



(14) 



where r\\ is the Caille parameter determined by the line 
shape. The experimental data for three different lipids 
indicated that Pfi could be represented by an exponen- 
tial exp(—a/Xfi), in agreement with the result of pertur- 
bation theory ||, but that A/; was significantly greater 
than 2A instead of exactly 2A given by perturbation the- 
ory. Since neither the perturbation theory nor the har- 
monic interpretation of the data are necessarily correct, 
it is valuable to test these predictions using simulations. 

Figure [l^ shows two ways of obtaining Pfi from the 
simulations. The first way uses the definition 



P = P 



fl 



Pb, 



(15) 



where P is the total osmotic pressure and Pb is the pres- 
sure with no fluctuations, i.e. for the membrane ex- 
actly in the middle of the space between the two walls 
with u(x,y) = 0. The second way uses Eq.([l3]). Fig .[To| 
shows that the simulated Pfi can be reasonably repre- 
sented by an exponential using either method of compu- 
tation, thereby supporting both theory and experiment. 
Either method gives decay lengths Xfi that exceed 2A, 
thereby supporting experiment. The two results for Pfi 
in Fig|l(] do not, however, agree perfectly, and the dis- 
crepancy grows for larger values of a. This is not sur- 
prising because the harmonic approximation is better for 
small a and progressively breaks down, especially when 
the bare potential no longer has a minimum at z = 0. 
This discrepancy suggests that one should expect some 
error when subtracting Pfi obtained from Eq.(|l3|) from 
P in Eq.(|l5|) to obtain Pi,, although the error is encour- 
agingly small. Nevertheless, future work in this direction 
can employ simulations to correct this discrepancy and 
to allow a better estimate of Pb from which Ph , A and H 
are obtained . 




a [A] 

FIG. 10. Simulation results for Pfi vs. a for A = 1, H = 4, 
K c = 0. 5 Il8j , A = 1.8A. Solid circles show Pfi obtained 
from Eq. (113) with a slope A/; = 4.1 A. Open circle show Pfi 
obtained from Eq.(|l3|) with a slope Xfi = 4.6A. 



VIII. CONCLUSIONS 

This paper solves accurately a model of constrained 
single membrane fluctuations. The new FMC simulation 
method provides a way to simulate accurately, with mod- 
est computer resources, the pressure and mean square 
fluctuation of a simple membrane between two hard walls 
with realistic potentials. This method is clearly superior 
to the more conventional PMC simulation method. Used 
with typical values of interaction parameters, it supports 
the idea of the exponential decay of fluctuational pres- 
sure, lending credibility to a simplified interpretation of 
X-ray scattering data in Q. Finally, the method, with 
minor modification, may be applied to studies of more 
complicated models, such as a stack of membranes or 
models of charged lipids and more sophisticated data 
analysis. 
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